home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / luc.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  8.5 KB  |  335 lines
  1. /* linalg/luc.c
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <string.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_vector.h>
  25. #include <gsl/gsl_matrix.h>
  26. #include <gsl/gsl_complex.h>
  27. #include <gsl/gsl_complex_math.h>
  28. #include <gsl/gsl_permute_vector.h>
  29. #include <gsl/gsl_blas.h>
  30. #include <gsl/gsl_complex_math.h>
  31.  
  32. #include <gsl/gsl_linalg.h>
  33.  
  34. /* Factorise a general N x N complex matrix A into,
  35.  *
  36.  *   P A = L U
  37.  *
  38.  * where P is a permutation matrix, L is unit lower triangular and U
  39.  * is upper triangular.
  40.  *
  41.  * L is stored in the strict lower triangular part of the input
  42.  * matrix. The diagonal elements of L are unity and are not stored.
  43.  *
  44.  * U is stored in the diagonal and upper triangular part of the
  45.  * input matrix.  
  46.  * 
  47.  * P is stored in the permutation p. Column j of P is column k of the
  48.  * identity matrix, where k = permutation->data[j]
  49.  *
  50.  * signum gives the sign of the permutation, (-1)^n, where n is the
  51.  * number of interchanges in the permutation. 
  52.  *
  53.  * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
  54.  * Elimination with Partial Pivoting).
  55.  */
  56.  
  57. int
  58. gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
  59. {
  60.   if (A->size1 != A->size2)
  61.     {
  62.       GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
  63.     }
  64.   else if (p->size != A->size1)
  65.     {
  66.       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
  67.     }
  68.   else
  69.     {
  70.       const size_t N = A->size1;
  71.       size_t i, j, k;
  72.  
  73.       *signum = 1;
  74.       gsl_permutation_init (p);
  75.  
  76.       for (j = 0; j < N - 1; j++)
  77.     {
  78.       /* Find maximum in the j-th column */
  79.  
  80.       gsl_complex ajj = gsl_matrix_complex_get (A, j, j);
  81.           double max = gsl_complex_abs (ajj);
  82.       size_t i_pivot = j;
  83.  
  84.       for (i = j + 1; i < N; i++)
  85.         {
  86.           gsl_complex aij = gsl_matrix_complex_get (A, i, j);
  87.               double ai = gsl_complex_abs (aij);
  88.  
  89.           if (ai > max)
  90.         {
  91.           max = ai;
  92.           i_pivot = i;
  93.         }
  94.         }
  95.  
  96.       if (i_pivot != j)
  97.         {
  98.           gsl_matrix_complex_swap_rows (A, j, i_pivot);
  99.           gsl_permutation_swap (p, j, i_pivot);
  100.           *signum = -(*signum);
  101.         }
  102.  
  103.       ajj = gsl_matrix_complex_get (A, j, j);
  104.  
  105.       if (!(GSL_REAL(ajj) == 0.0 && GSL_IMAG(ajj) == 0.0))
  106.         {
  107.           for (i = j + 1; i < N; i++)
  108.         {
  109.           gsl_complex aij_orig = gsl_matrix_complex_get (A, i, j);
  110.                   gsl_complex aij = gsl_complex_div (aij_orig, ajj);
  111.           gsl_matrix_complex_set (A, i, j, aij);
  112.  
  113.           for (k = j + 1; k < N; k++)
  114.             {
  115.               gsl_complex aik = gsl_matrix_complex_get (A, i, k);
  116.               gsl_complex ajk = gsl_matrix_complex_get (A, j, k);
  117.                       
  118.                       /* aik = aik - aij * ajk */
  119.  
  120.                       gsl_complex aijajk = gsl_complex_mul (aij, ajk);
  121.                       gsl_complex aik_new = gsl_complex_sub (aik, aijajk);
  122.  
  123.               gsl_matrix_complex_set (A, i, k, aik_new);
  124.             }
  125.         }
  126.         }
  127.     }
  128.       
  129.       return GSL_SUCCESS;
  130.     }
  131. }
  132.  
  133. int
  134. gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
  135. {
  136.   if (LU->size1 != LU->size2)
  137.     {
  138.       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
  139.     }
  140.   else if (LU->size1 != p->size)
  141.     {
  142.       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
  143.     }
  144.   else if (LU->size1 != b->size)
  145.     {
  146.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  147.     }
  148.   else if (LU->size2 != x->size)
  149.     {
  150.       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
  151.     }
  152.   else
  153.     {
  154.       /* Copy x <- b */
  155.  
  156.       gsl_vector_complex_memcpy (x, b);
  157.  
  158.       /* Solve for x */
  159.  
  160.       gsl_linalg_complex_LU_svx (LU, p, x);
  161.  
  162.       return GSL_SUCCESS;
  163.     }
  164. }
  165.  
  166.  
  167. int
  168. gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
  169. {
  170.   if (LU->size1 != LU->size2)
  171.     {
  172.       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
  173.     }
  174.   else if (LU->size1 != p->size)
  175.     {
  176.       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
  177.     }
  178.   else if (LU->size1 != x->size)
  179.     {
  180.       GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
  181.     }
  182.   else
  183.     {
  184.       /* Apply permutation to RHS */
  185.  
  186.       gsl_permute_vector_complex (p, x);
  187.  
  188.       /* Solve for c using forward-substitution, L c = P b */
  189.  
  190.       gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);
  191.  
  192.       /* Perform back-substitution, U x = c */
  193.  
  194.       gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);
  195.  
  196.       return GSL_SUCCESS;
  197.     }
  198. }
  199.  
  200.  
  201. int
  202. gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual)
  203. {
  204.   if (A->size1 != A->size2)
  205.     {
  206.       GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
  207.     }
  208.   if (LU->size1 != LU->size2)
  209.     {
  210.       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
  211.     }
  212.   else if (A->size1 != LU->size2)
  213.     {
  214.       GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
  215.     }
  216.   else if (LU->size1 != p->size)
  217.     {
  218.       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
  219.     }
  220.   else if (LU->size1 != b->size)
  221.     {
  222.       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
  223.     }
  224.   else if (LU->size1 != x->size)
  225.     {
  226.       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
  227.     }
  228.   else
  229.     {
  230.       /* Compute residual, residual = (A * x  - b) */
  231.  
  232.       gsl_vector_complex_memcpy (residual, b);
  233.  
  234.       {
  235.         gsl_complex one = GSL_COMPLEX_ONE;
  236.         gsl_complex negone = GSL_COMPLEX_NEGONE;
  237.         gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, residual);
  238.       }
  239.  
  240.       /* Find correction, delta = - (A^-1) * residual, and apply it */
  241.  
  242.       gsl_linalg_complex_LU_svx (LU, p, residual);
  243.  
  244.       {
  245.         gsl_complex negone= GSL_COMPLEX_NEGONE;
  246.         gsl_blas_zaxpy (negone, residual, x);
  247.       }
  248.  
  249.       return GSL_SUCCESS;
  250.     }
  251. }
  252.  
  253. int
  254. gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
  255. {
  256.   size_t i, n = LU->size1;
  257.  
  258.   int status = GSL_SUCCESS;
  259.  
  260.   gsl_matrix_complex_set_identity (inverse);
  261.  
  262.   for (i = 0; i < n; i++)
  263.     {
  264.       gsl_vector_complex_view c = gsl_matrix_complex_column (inverse, i);
  265.       int status_i = gsl_linalg_complex_LU_svx (LU, p, &(c.vector));
  266.  
  267.       if (status_i)
  268.     status = status_i;
  269.     }
  270.  
  271.   return status;
  272. }
  273.  
  274. gsl_complex
  275. gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
  276. {
  277.   size_t i, n = LU->size1;
  278.  
  279.   gsl_complex det = gsl_complex_rect((double) signum, 0.0);
  280.  
  281.   for (i = 0; i < n; i++)
  282.     {
  283.       gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
  284.       det = gsl_complex_mul (det, zi);
  285.     }
  286.  
  287.   return det;
  288. }
  289.  
  290.  
  291. double
  292. gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
  293. {
  294.   size_t i, n = LU->size1;
  295.  
  296.   double lndet = 0.0;
  297.  
  298.   for (i = 0; i < n; i++)
  299.     {
  300.       gsl_complex z = gsl_matrix_complex_get (LU, i, i);
  301.       lndet += log (gsl_complex_abs (z));
  302.     }
  303.  
  304.   return lndet;
  305. }
  306.  
  307.  
  308. gsl_complex
  309. gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
  310. {
  311.   size_t i, n = LU->size1;
  312.  
  313.   gsl_complex phase = gsl_complex_rect((double) signum, 0.0);
  314.  
  315.   for (i = 0; i < n; i++)
  316.     {
  317.       gsl_complex z = gsl_matrix_complex_get (LU, i, i);
  318.       
  319.       double r = gsl_complex_abs(z);
  320.  
  321.       if (r == 0)
  322.     {
  323.       phase = gsl_complex_rect(0.0, 0.0);
  324.       break;
  325.     }
  326.       else
  327.         {
  328.           z = gsl_complex_div_real(z, r);
  329.           phase = gsl_complex_mul(phase, z);
  330.         }
  331.     }
  332.  
  333.   return phase;
  334. }
  335.